program chem_reaction_test

   use mod_chem
   use stdlib_string_type
   use iso_fortran_env, only: dp => real64
   implicit none

   type(type_chem_phases) :: phases
   integer, parameter :: n = 1981
   real(dp), allocatable :: T(:), X(:, :), rho_compared(:), H, U
   real(dp) :: p = 101325.0_dp, rho
   real(dp), allocatable :: molar_density(:), creation_rate(:), destruction_rate(:), net_rate(:)

   integer :: l, i, ll

   integer :: file_io(100)

   phases = chem_file_load("air_5species_NASA9.yaml", "air")
   allocate (T(n), X(5, n), rho_compared(n))
   open (newunit=file_io(100), file="5S_P_xi.plt")
   read (file_io(100), *)
   read (file_io(100), *)
   do l = 1, n
      read (file_io(100), *) T(l), X(2, l), X(5, l), X(1, l), X(3, l), X(4, l), rho_compared(l)
   end do
   close (file_io(100))

   l = 1981

   call phases%set_state(T=T(l), X=X(1:5, l), P=p)
   print *, phases%get_mean_molecular_weight()
   do i = 1, phases%n_reactions
      write (*, '(1X, 2A)') phases%reactions(i)%equation, ":"
      write (*, '(1X, 2A)') "type: ", phases%reactions(i)%type
      write (*, '(1X, 3(A, G20.7))') "T:", phases%get_temperature(), "P:", phases%get_pressure(), "D:", phases%get_density()
      write (*, '(1X, A, 5G20.7)') "molar_density:", phases%get_molar_density()
      molar_density = phases%get_molar_density()
      write (*, '(1X, A, 5G20.7)') "fwd elementary coeffs:", &
         (molar_density(ll)**phases%reactions(i)%sto_num_fwd(ll), ll=1, phases%n_species)
      write (*, '(1X, A, 5G20.7)') "rev elementary coeffs:", &
         (molar_density(ll)**phases%reactions(i)%sto_num_rev(ll), ll=1, phases%n_species)
      if (phases%reactions(i)%type == "three-body") then
         write (*, '(1X, A, 2G20.7)') "three-body coeff", sum(molar_density*phases%reactions(i)%efficiencies)
         write (*, '(1X, A, 5G20.7)') "effiecients:", phases%reactions(i)%efficiencies
      end if
      write (*, '(1X, A, G20.7)') "k:", phases%reactions(i)%Arrhenius(T(l))
      write (*, '(1X, A, G20.7)') "kf:", phases%reactions(i)%get_fwdrate_constant(T(l))
      write (*, '(1X, A, G20.7)') "equilibrium constant kp: ", phases%reactions(i)%get_equilibrium_constant_p(T(l))
      write (*, '(1X, A, G20.7)') "equilibrium constant kx: ", phases%reactions(i)%get_equilibrium_constant_x(T(l))
      write (*, '(1X, A, G20.7)') "kb:", phases%reactions(i)%get_revrate_constant(T(l))
      write (*, '(1X, A, G20.7)') "fwd rate: ", phases%reactions(i)%get_fwdrate_progress(phases%get_molar_density(), T(l))
      write (*, '(1X, A, G20.7)') "rev rate: ", phases%reactions(i)%get_revrate_progress(phases%get_molar_density(), T(l))
      write (*, '(1X, A, G20.7)') "net rate: ", phases%reactions(i)%get_netrate_progress(phases%get_molar_density(), T(l))
      write (*, '(1X, A, 5G20.7)') "creation rate", &
         phases%reactions(i)%get_creation_production_rates(phases%get_molar_density(), T(l))
      write (*, '(1X, A, 5G20.7)') "destruction rate", &
         phases%reactions(i)%get_desctruction_production_rates(phases%get_molar_density(), T(l))
      write (*, '(1X, A, 5G20.7)') "net rate", phases%reactions(i)%get_net_production_rates(phases%get_molar_density(), T(l))
   end do
   creation_rate = phases%get_creation_production_rates()
   destruction_rate = phases%get_destruction_production_rates()
   net_rate = phases%get_net_production_rates()
   write (*, *)
   write (*, *) "creation/destruction/net rate of productions:"
   do i = 1, phases%n_species
      write (*, '(1X, A, 3G20.7, A)') phases%species(i)%name, creation_rate(i), destruction_rate(i), net_rate(i), "   kmol/m^3/s"
   end do

   open (newunit=file_io(1), file="equ_constant.plt")
   open (newunit=file_io(2), file="fwd_rev_rate_constant.plt")
   open (newunit=file_io(3), file="fwd_rev_rate.plt")
   open (newunit=file_io(4), file="fwd_rev_rate_new.plt")
   open (newunit=file_io(5), file="production_rate.plt")
   do l = 1, 5
      write (file_io(l), '(A)') "variables = T"
   end do
   do i = 1, phases%n_reactions
      write (file_io(1), "(1X, A)") "Eq"//char(string_type(i))//"_Equ_c"
      write (file_io(2), "(1X, 3A)") "Eq"//char(string_type(i))//"_fwd_rate_constant", new_line(''), &
         "Eq"//char(string_type(i))//"_rev_rate_constant"
      write (file_io(3), "(1X, 3A)") "Eq"//char(string_type(i))//"_fwd_rate", new_line(''), &
         "Eq"//char(string_type(i))//"_rev_rate", new_line(""), &
         "Eq"//char(string_type(i))//"_net_rate"
   end do

   do i = 1, phases%n_species
      write (file_io(4), "(1X, A)") "Eq"//char(string_type(i))//"_fwd_rate"
   end do
   do i = 1, phases%n_species
      write (file_io(4), "(1X, A)") "Eq"//char(string_type(i))//"_rev_rate"
   end do
   do i = 1, phases%n_species
      write (file_io(4), "(1X, A)") "Eq"//char(string_type(i))//"_net_rate"
   end do

   write (file_io(5), "(1X, A)") (phases%species(i)%name//"_creation_rate", i=1, phases%n_species), &
      (phases%species(i)%name//"_destruction_rate", i=1, phases%n_species), &
      (phases%species(i)%name//"_net_rate", i=1, phases%n_species), &
      (phases%species(i)%name//"_mass_creation_rate", i=1, phases%n_species), &
      (phases%species(i)%name//"_mass_destruction_rate", i=1, phases%n_species), &
      (phases%species(i)%name//"_mass_net_rate", i=1, phases%n_species)

   do l = 1, n
      call phases%set_state(T=T(l), X=X(1:5, l), P=p)
      write (file_io(1), '(1X, 8G16.8)') T(l), (phases%reactions(i)%get_equilibrium_constant_x(T(l)), i=1, phases%n_reactions)
      write (file_io(2), '(1X, 20G16.8)') T(l), (phases%reactions(i)%get_fwdrate_constant(T(l)), &
                                                 phases%reactions(i)%get_revrate_constant(T(l)), i=1, phases%n_reactions)
      write (file_io(3), '(1X, 20G16.8)') T(l), &
         (phases%reactions(i)%get_fwdrate_progress(phases%get_molar_density(), T(l)), &
          phases%reactions(i)%get_revrate_progress(phases%get_molar_density(), T(l)), &
          phases%reactions(i)%get_netrate_progress(phases%get_molar_density(), T(l)), i=1, phases%n_reactions)
      write (file_io(4), '(1X, 20G16.8)') T(l), &
         (phases%reactions(i)%get_fwdrate_progress(phases%get_molar_density(), T(l)), i=1, phases%n_reactions), &
         (phases%reactions(i)%get_revrate_progress(phases%get_molar_density(), T(l)), i=1, phases%n_reactions), &
         (phases%reactions(i)%get_netrate_progress(phases%get_molar_density(), T(l)), i=1, phases%n_reactions)
      write (file_io(5), '(1X, 40G16.8)') T(l), phases%get_creation_production_rates(), &
         phases%get_destruction_production_rates(), phases%get_net_production_rates(), &
         phases%get_creation_production_rates()*phases%get_molecular_weights(), &
         phases%get_destruction_production_rates()*phases%get_molecular_weights(), &
         phases%get_net_production_rates()*phases%get_molecular_weights()
   end do

   do l = 1, 5
      close (file_io(l))
   end do

end
